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A a eee nl one 


Abstract 


it is shown that, for nuclear physical reasons, an end 
velocity higher than 0.134c cannot be attained by a rocket 
ising a propulsion system which maximizes the end velocity, 
vt given energy content and end mass, by varying the ex- 
aust velocity or by using propellants with variable com- 
dustion heat. 

Limiting Factors for the End Velocity of a 
Rocket 


In discussions about reaching high velocities with a 
“ocket and especially reaching relativistic velocities 
1}-{6] [10] [12], two limiting factors are outstandingly 
spparent. 

‘One of these is that the mass conversion factor gt 

# nuclear fuel, (that fraction of fuel mass which can be 
sonverted into energy) is low. As a consequence un- 
ealistic values of initial mass to final mass would be 
‘equired to obtain such high end velocities. 
The other is that the specific power of the rocket 
mgine and the burning time are limited. This second 
actor can be expressed in the form: the energy F for 
: given end mass /, of the rocket, or the product 


E/M. = P/M,-tis limited, 


vhere P/M, is the power engine mass and ¢ the burning 
ime. (Here we equalize the end mass and the engine 
nass, that means we include the payload and structure 
mass in the engine mass for simplicity.) 

The end velocity which cannot be surpassed is 
letermined by H/M, or by this product (P/M.)t as 
nay be seen from the formula 

Mo cies 
2 
vhere wu, is the end velocity. The “characteristic velocity 
= 4/2E/M.” obtained in transforming the whole 
nergy into kinetic energy of the end mass cannot be 
urpassed. 

To give an example of how our present technology 
ompares with the task of reaching relativistic veloci- 
ies, it may be mentioned that a burning time of 50 
ears would be required, using an optimal propulsion 
ystem and an engine having the specific power of the 


* Adapted from Colloquim notes 20 Jan. 1958, AFMDC, 
folloman A. F. B. 

7 N.A.S.A. Langley Research Center. 

tg has no relationship to the acceleration in the gravita- 
ional field of the earth for which normally the same letter is 
sed. 
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V2 (here assumed to be 500 kw/kg engine mass) in 
order to attain an end velocity of only 0.134¢c. Such 
an optimum propulsion system is defined in the fol- 
lowing discussion. 


An Optimum Propulsion System with Variable 
Exhaust Velocity 


Propulsion systems with constant exhaust velocity 
Cp* (Cp = Cpsopellant, relative to the vehicle) are not the 
most favorable systems with respect to end velocity, 
if energy # and end mass M, are given. An optimum 
propulsion system can be visualized, where the exhaust 
velocity c, as a function of the decreasing mass of the 
rocket is varied in such a way that the end velocity u, 
of the rocket will be a maximum [7]-[11]. Using classical 
mechanics the exhaust velocity must be inversely 
proportional to the instantaneous mass of the rocket. 
The relativistic treatment of the problem was given by 
H. Knothe [10] using a general inequality, which con- 
tains Schwartz’s inequality as a special case. The above 
classical approximation for the dependence of ¢c, on 
the rocket mass is a good approximation at the rela- 
tivistically low velocities considered here. We may sum- 
marize some of the qualities of this propulsion system 
[6b]. 

The end velocity u. which can be approached by 
increasing the mass ratio W,/M, (M, initial mass of 
the rocket) is in classical approximation u. = »/2E/M,; 
this means there is nearly no loss of energy in the ex- 
hausted jet. A predetermined end velocity u, is reached 
with a minimum of energy, that is LE = M,(u,’/2), 
if M,/M,. > 1. Thus the optimum system has the 
advantage that the product of specific power of the 
engine and burning time for attaining the _ prede- 
termined end velocity will be a minimum since this 
product (P/M,)t is proportional to the total energy 
E = P-t. Assuming a constant power output of the 
propulsion unit the acceleration will be constant; 
however this is not true of the thrust. Thus the opti- 
mum system avoids the high acceleration of the c, = 
constant system at the end of the propulsion phase. 


Maximal Final Velocity of the Optimum System in 
Comparison to that of the System having Constant 


* The reader is reminded that c, is the exhaust velocity and 
is not the specific heat at constant pressure which usually 
carries the same symbol. 
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Exhaust Velocity if Specific Power « Burning Time 

is Given 

For quantitative comparison with a propulsion system 
having constant exhaust velocity the classical ap- 
proximations for u. may be recalled [10]. If cp is any 
function of the variable mass M of the rocket, the 
Schwartz’s inequality can be applied for comparison of 


Ae a7 aM (1) 
and 
2 
B= / = aM (2) 
as follows: 


ue =({5 aM) < [ee am { 7,aM 


The equality (Optimum System) holds only if 
Cp = \/M 
where ) is related to H and M,, M, by 


ia pes Gey Bi) 1 
ae tau = (2-4), 


2H 

bes 
JE eas 
The combustion heats H(J/) of the variable fuel have 
to be 


(3) 
where the average combustion heat H is defined by 


E = H(M, — M.). Hence from (1) it follows that for 
the optimum system (c, ~ variable): 


p. ed ae Tey a eae 
Ute i Toe 4/ 28 ap ir) 

= (M, neu. 
4/ a r i\(1 ir im) ~ 
M, M, 


Considering the c, = const. system we have the 
known equations: 


I 
5 
La 


2 6 
H =c, /2 = const. 


Uy = I65 In; e = \/2H 


= H(M, — M,) 


ie 


we Pedy ise pee 
M./M, — ie 
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(B) 


If the engine and the mission is power- and tin 
limited, that means if 


E Ds 
M. M. 


then a maximum increase of 25 percent for uw. can 
obtained by using the optimum system, as is fi 
shown by D. B. Langmuir [8] and T. H. Irving | 
and can be recognized from (A) and (B) as folloy 

For the cp = const. system, we see from (B), varyi 
the mass ratio, that at a given characteristic veloc 

= »/2E/M, the end velocity ue has a maximt 
value u. = 0.8-v- at the mass ratio M,/M,. = 5 a 
at the exhaust velocity cp = 0.5-v.. 

For the optimum system (cp ® 1/M), u. in ( 
approaches its maximum u. = ~/2E/M, as M,/. 
is increased. Practically the maximum, the chi 
acteristic velocity itself, is reached at M,/M, = 
For a given characteristic velocity v. = 0.134¢ be 
the systems are compared in Fig. 1 (dashed line 


it is giv 


Comparison of the Optimum Propulsion Sy 
tems with the System having Consta 
Exhaust Velocity, if a Low Mass Co 
version Factor of the Fuel Matter 
Given 


Let us assume, however, that the engine and t 
mission are practically not time limited or that t 
specific power of advanced engines for long burni 
times could be essentially increased, so that the prodt 
specific power x burning time yields a characteris 
velocity essentially larger than 0.134 e. 

Then the question arises of how the end veloci 
of a rocket using such an optimum propulsion syst 
will be limited by the low mass conversion factor 
nuclear fuel, that is to say which end velocities can 
attained with the c, ~ 1/M system under the ec 
dition that the energy unit is coupled with a gre 
amount of mass. Most of the work considering re 
tivistic mechanics and propulsion has dealt with t 
Cp = const. system [1]-[6al [12]. The relativistic fi 
mulas for u. for such systems are given by Ackeret | 
Saenger [4] ete. [2]-{6] [10] [12]. From these formu 
can be calculated the mass conversion factors whi 
would be needed to reach relativistic velocities, 
several given mass ratios. Referring to the mass cc 
version factor, it is emphasized in [6a] that there exi 
an upper limit for the velocity attainable with a c, 
const. system, at a few tenths of the velocity of lig 
using normal matter as fuel. Nature put a basie barr 
for the conversion of matter into energy, this barr 
is the stability of protons, or the law of conservati 
of the number of protons and neutrons during a 
nuclear process which involves normal matter at te 
peratures less than 10“°K. Only the mass equival 
of the binding energy of the nuclei is available 
conversion into heat or kinetie energy in the nucl 
processes of fission and fusion and the mass convers’ 


actor for any such process cannot surpass the value of 
|= 1000; a8 is shown in the experimentally established 
at ding - Energy - Per - Particle curve (BEPP - Curve). 
ere IS seen, at present, no possibility of creating 
aacroscopic amounts of anti-matter on earth and to 
ave a fuel with a higher mass conversion factor then 
— roo0, disregarding the problems of controlling the 
iatter and anti-matter reaction and direction of the 
nnihilation radiations. 
The upper limit of the end velocity of a rocket 
sing the described optimum propulsion system may 
e calculated, it no fuel with a mass conversion factor 
Jarger than +259 is available.* 
Thus we have a problem, which differs from that 
reated in the above section where H/M, was given. 
Ve compare the c, = woes system with the c, ~ 1/M 
ystem if @max = Hmax/c’ is given as the same for both 
ystems, where H is the combustion heat of the fuels. 
The optimum propulsion system requires a fuel with 
‘ariable heat of combustion and this variable range 
nust increase as the assumed mass ratio of the rocket 
ncreases. In the classical approach the limits for the 
equired mass conversion factors g of the fuel are 
equation 3): 


vii ae eee 
ene eee oar 


Jmax = = 


rhere H, the average heat of combustion, is defined by 


= 5 E 
et 4 
aT SM, 4) 
lowever, the maximum mass conversion factor 
a 
bit 9 (5) 


eee aT 1000 


f normal matter, which is the only possible fuel, 
annot be surpassed. Therefore, the average heat of 
ombustion of the fuel will decrease as the assumed 
aass ratio increases. The fuel degenerates to a certain 
egree, although the attainable end velocity increases. 
‘he energy content of such rockets is limited by the 
alue ¢’:max:M, for large M./M., as can be seen by 
ubstituting gmax(M./M.) for g from equation (5) in 
quation (4). We obtain in equation (4) 


3 2 M, 
a H(M, =) M.) = Gnas i, (1 — | 
arn) 
& Jmax'¢ -M, for large: ir. 
his energy content cannot be increased by increasing 


he mass ratio or the mass of the fuel as in the case of 


* We consider here only the case of variable c, treated up to 
ow, that is the case where the propellant is the combustion 
roducts of the fuels or where the combustion products having 
‘st mass of a separate energy source are contained in the 
ropellant. The system with variable c,, where the combustion 
roducts of a separate energy source are partially thrown 
yerboard without contributing to the thrust, to reach higher 
chaust velocities, is treated in a following note. 


propulsion with cp, = const. Therefore, supposing 


9 : 
Imax = 7000’ the greatest end velocity, that can be 


approached is given by the formula 
= V/2E/M,. = V2CGmax = 0.134 
which results by inserting # from (6), see also equation 
(A). 
For comparing quantitatively the c, = const and 


Cp» ~ 1/M systems at given gmax, We consider again 
the equations (A) (B) for u., introducing the maximum 


mass conversion factor gmax = or the maximum 


9 
1000 
combustion heat Himax = Jimax*C- 

For the c, ~ 1/M system we get from (A) and 


H = Jmax:¢ (M./M.) 
at M, M,. 
Ue = /9 Ee ee: 
V 2H 4/ GE we iz) 


or Ss M, 
= N/ 2 ies (1 a ae 


—>~/2H max for increasing M,/M,. 


For cy = const however, we obtain from (B) 
= / 2H In Mi/M,, 


which increases arbitrarily if 17,/M, increases. 

We have thus the surprising result that if gmax is 
given the c, = const propulsion system is far superior 
at large mass ratios in respect to reaching high end 
velocities. The maximum mass conversion factor 
g = te restricts the applicability of the optimum 
method to an end velocity smaller then 0.134c. Using 
the optimum propulsion at the beginning, and then 
continuing the propulsion with gmax = const, a higher 
mass ratio would be required than with c, = const 
propulsion using the best fuel available from the 
beginning on. In the latter case of course a greater 
amount of the optimum fuel and of energy would be 
consumed. In all these considerations we assume that 
the internal efficiency of the propulsion motor is equal 
to one, and furthermore, that such an optimum thermo 


could be realized and 


9 
1000 
controlled in a rocket. This problem is yet unsolved 
even for the fusion reaction D + T — H,4 + n having 


nuclear reaction with g = 


rang Using the rela- 
tivistic kinematic equations only minor corrections 
will result at these relativistically low velocities of the 
rocket and propellant. 

In Fig. 1 the solid lines show the values u,/c, which 
can be attained by varying the mass ratio for both the 


the mass conversion factor g = 


systems: 


Cp ™ var 
M if Omens = i ae 00 is given. 


Cp = const | 


= 
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Fig. 1. A comparison of the optimum system (cp variable) 
and the constant c, system with reference to attainable end 
velocities. 


— For a given maximum mass conversion factor, 


9 


Bese 600 


: E Power 
—-—-— For a given = 


- X burning time = Hnax 
M,. engine mass 


= Emax X ¢? which is a given characteristic velocity, 


LG Te 
M. 


e 


The shadowed area would be only accessible by cp = 
const. The dotted area cannot be reached at all if g < 
100° Especially a genuine photon rocket would pre- 
suppose a mass conversion factor one. 

The dashed lines show the upper limits of end 
velocities which could be attained, if H/M, = 
(P/M.)t Jmax*C? = 8.1 X 10 erg/gr, which cor- 
responds to a characteristic velocity 0.134c. In this 
case, where the amount of energy for an end mass 
M, 1s given, the c» = const system is inferior. For the 

: an ¢ 
Cp ~ variable system the conditions gmax = den and 


1000 
E/M, = 8.1 X 108 = “gnax”’-C are identical for large 
M,/M., of course for a low mass ratio the former 
condition is more restrictive. 
It should, however, be emphasized that at given 
9 Siem : : ; 
2max = 1000 rhe » miss and engine time is not 


limited—the ce, = constant system delivers an es- 
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sentially higher u. and does not imply any restrictio 
the available energy or end velocity except that of 
quiring that the initial-to-end mass ratio be feasi 
We can thus compensate to a certain degree for 


. . 9 
relatively low exhaust velocity given by g = 1000 


using a high mass ratio. Of course, however, we can 
hope to attain relativistic velocities resulting in 
essential time dilatation having such low mass 4 

2 


version factor of the fuel. 
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Attitude Drift of Space Vehicles 


4 William 'T. Thomson* and Gordon S. Reitert 


Abstract 


_ The attitude of a body of revolution spinning in the absence 
ee forces is not a constant when energy dissipation 
takes place. Elastic vibration, induced by gyroscopic action, 
results in a dissipation of energy and a change in the pre- 
cession cone angle @. This paper examines the effect of energy 
dissipation on the spinning body and evaluates the time 
required for a body of given configuration to undergo a 
specified change in attitude. 


I. Introduction 


_ Several proposals have been studied for changing the 
attitude of spinning space vehicles by using active 
internal control torques. [1] Attitude changes can also 
occur in a completely passive system without external 
moments, since the attitude of a body of revolution 
spinning in the absence of external forces is not a con- 
stant when energy dissipation takes place. This principle 
has been used in several devices which stabilize spin- 
ning satellites by intentional dissipation of energy. [2], 
3] 

Elastic vibration, induced by gyroscopic action, also 
results in a dissipation of energy and a change in the 
precession cone angle 6. This paper examines the effect 
of energy dissipation due to elastic vibration on a spin- 
ning body, and evaluates the time required for a body 
of given configuration to undergo a specified change in 
uttitude. 

The moment-free motion of an unsymmetric body 
with principal moments of inertia A, B, C, is an un- 
steady periodic precession and nutation about the 
esultant angular momentum vector / fixed in space. 
Steady rotation is possible only about the principal 
ixis of maximum or minimum moment of inertia, the 
wincipal axis of intermediate moment of inertia being 
mstable. 

For a body of revolution A, A, C, the moment-free 
notion is a steady precession of the spin axis at a 
constant angle 6, about the resultant angular momen- 
uum vector / fixed in space. Steady rotation is again 
Jossible about the axis of maximum or minimum 
noment of inertia, and the principal axis of intermedi- 
ite moment of inertia does not exist. 

In either case, the axis of maximum or minimum 
noment of inertia is considered to be stable in that if 
he spin axis deviates slightly from the resultant an- 

* Professor of Engineering, University of California, Los 
ngeles. Consultant, Space Technology Laboratories. 

+ Member of the Technical Staff, Space Technology 
aboratories, Los Angeles 45, Calif. 


gular moment vector, there is no tendency of this 
deviation to grow. This statement is true only for a 
perfectly rigid body in the absence of external moment. 

In an elastic body, deformation between particles 
will always take place, resulting in some dissipation of 
energy. When the dissipation of energy is taken into 
account, we must revise our statement of stability in 
that it is possible for a small deviation of the spin axis 
to grow into a large one and eventually result in a 
complete change in attitude of the body. For such 
bodies, only the principal axis of maximum moment of 
inertia is stable, and the axis of minimum moment of 
inertia is one of unstable equilibrium. 

These facts were actually observed in the Explorer I 
Satellite, [4] which was spin stabilized about the longi- 
tudinal axis of minimum moment of inertia. The 
flexible antennas of the satellite provided an excellent 
source for energy dissipation, and in one revolution 
around its orbit (approximately 90 minutes) the 
Explorer I was observed to be tumbling at an attitude 
of 6 = 60° instead of spinning about its longitudinal 
axis at 6 = 0. The remedy for this behavior is obviously 
to shorten the longitudinal dimensions of the satellite 
so that the moment of inertia about the longitudinal 
spin axis is greater than that about the transverse 
pitch or yaw axis. However, the problem still exists for 
missiles which are long slender bodies and inherently 
unstable about the spin axis. Here the important 
question is how long can the spinning missile coast in 
a moment-free condition without an appreciable change 
in its attitude. 


II. Energy Considerations of Stability 


We will now examine the basis for stability from an 
energy point of view. For a body of revolution with 
principal moments of inertia A, A, C, as shown in 
Figure 1, the moment-free motion is that of steady 


precession described by the equations 
CIES ee oes 

(A —C)cos@ Acosé (1) 
ws = 6 + ¥ cos 8 


v 


Since the moment is zero the angular momentum vec- 
tor fh is a constant with the following components 
along the body axes 


h=A (a3 -F wo] ) + Cask (2) 
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Fig. 1. Coordinate System of Body Axes 1, 2, 3 


The magnitude of h can be specified in terms of its 
initial attitude 6& 0, and assuming 6 to be small, 
the following relationships must hold. 


h — Cao 
@3 = Wo COS 6 
; 6) 
y ae A Wo (3) 


I 


d (1 — ae cos 0 


We must examine the kinetic energy of rotation 
T= 3A (Gir ale we’) + 1Ows" 
For 6 small, we have 
AVwor + wf = hsin 6 
Caw; = h cos 0 


(4) 


so that 7’ can be expressed in the form 


Ig it ar C re 
T = 5 Cw E + (4 - 1) sin” | (5) 


Equation (5) indicates that @ remains constant 
provided T is constant. However with dissipation of 
energy 7’ must decrease. Differentiating (5) we have, 
as in References [5] and [6] 


Y 
ri 2 fC ; ; 
10 = (Ola (S - 1) sin 0 60804 (6) 
and since 7’ must always be negative, 6 is negative for 


(C/A) > 1 and positive for (C/A) < 1. Thus the 
principal axis of minimum moment of inertia is one of 
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unstable equilibrium, and a small deviation of t 
spin axis will increase due to energy dissipation whi 


(C/A). <1. 


III. Dissipation of Energy 


Assuming an elastic body, the energy dissipated p 
unit volume per cycle of stress can be assumed to be 


: 


ve 

2E ore | 
where y is a hysteretic damping factor establishing tl 
fraction of the elastic energy which is dissipated : 
shown by the shaded area in Fig. 2. Dividing by tl 
time ¢ per cycle of stress, and integrating over tl 
entire structure, the rate of energy dissipation can |] 
found. Thus the equation to be solved is of the gener, 
form 


Y 


; 
es / oEe dV = Cao & — 1) sin @cos66 (& 

In examining the source of cyclic stressing, fre 
vibration can be discarded since it will soon damp ou 
Steady cycling of stress is however induced by tl 
gyroscopic precession, and these stresses are repeate 
at the rate ¢ and 2¢, as we will presently show. 

The excitation for the cyclic stress is the acceler: 
tion. Choosing an arbitrary point on the structure an 
orienting the plane 1, 0, 3, through it, the displaceme 
vector for the specified point is 


F= H+ 2k € 


With 6 small, the angular velocity and acceleration | 
the coordinate axes 1, 2, 3 are 


® = ysin Osing? + Ysin 6cos¢j + (6+ Weos yk (1 


Fig. 2. Energy Dissipated by Hysteretic Damping 


A @ Sone 
use. [ee +aoXo 
dt i,j,k,=const CLT) 


= ¢W sin 6 (cos 7 — sin $7) 


f= +a + HX (@XF)+aX Fa 20 XK 0 (12) 


the result after some algebraic reduction is 

a= |-£(¢ + W’) + W sin’ Osin’ ¢ 

Fe — 2&dy cos 6 + 2f” sin 0 cos @ sin |i 
? + (i sin’ 6 sin ¢ cos @ 

> + zy’ sin 6 cos 6 cos $) J 
~ + (2&¢y sin 6 sin d 

+ £7 sin 0 cos @ sin @ — 2’ sin’ O)k 


(13) 


4 somewhat more convenient form of (13) results by 
eliminating ¢ and ¥ by means of (3). 


7 y\ 2 
Gd = wy | - é +¢é (5) sin’6 sin’ 
A A 
2 y\ 2 
(9) — 1 eos 6+ 2 (“) sin @ cos @sin 6 


2 CPS his 1 (Ones 
+ wo | (| —] sin’ @sin dcos ¢ + z 7 sing (14) 


-COs 8COS 6) + w'| - ie _ 2) sin 0 
Ca Vata ao 
‘cos 8 sin @ — z a sin 6\k 


Since the only time varying quantity in (14) (assum- 
ng 6 to be negligible) is ¢ = ¢f, it is evident that the 
yclic stress is repeated at a rate ¢ and 2¢. It should be 
pointed out that for slender bodies like missiles, C/A 
$ small compared to unity and the predominant 
ariable acceleration term is 


Gi, = 2wy £ So sin # cos 6sin ¢: k (15) 


hich is repeated in the time 


2 py) 
—— = ot iwhen C <A 
1- zh wp COSO 
; (16) 
2Qr 


= when C>A 


¢ — 1) Wo cos 6 


As an example of the simplest kind, we will consider 
wo solid disks connected by a flexible tube as shown 


\V. Example 


Fig. 3. Satellite Configuration, Displacement and Mo- 
ment Distribution. 


in Fig. 3. We will let C; and A, be the moment of inertia 
of each disk about its own polar and diametric axes. 
The gyroscopic moment required by each disk is 


Mg = Ci(¢ + cos 6)¥ sin 6 — AW sin 6 cos 8 (17) 


Since the moments of inertia about the center of mass 
of the body are 


C= 264 
A > 2(Ay + ml’) 
Equation (17) can be rewritten as 


Me = 4(C(¢ + ¥ cos 0) sin 6 — Ay sin 6 cos 6] 
+ ml'V sin 0 cos 6 


The first term, however, is the moment about the 
center of mass which is zero and from which (1). is 
obtained. We are thus left with 


Ma, = mW sin 0 cos 0 
(18) 


Fl cos 0 
where F = my’ sin 6 is the centripetal force of the 
precessing disk. 

The effect of MW, and F on the flexible tube is shown 
in Fig. 3. At point z along the tube, measured from 
the center of mass, the bending moment is 


Meee Mey (19) 


and the expression for the maximum stress becomes 


M,& 
Vs 


ee =m Cy : 3 sin 6 cos 0 


it (eV ccat @ © kin B.0O8 0 
= ) m A Wo / (tes COs 


which is repeated at the rate given by (16). The rate 
of energy dissipation as given by the left side of (8) 


(20) 


is then 


2 2 y¥\ 4 Y e 
ms ss) V Gy G — 1) wy sin’@ cos 6 (21) 
TL ZL i 


and the rate of change of the attitude angle @ becomes 
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@ (DEGREES) 


Fig. 4. Variation in the Rate of Tumbling 


7 a V (i) ae ; 
== = 9 Sin 6 cos 6 


= K sin 6 cos 0 


A plot of (22) is shown in Fig. 4. Since 6 is zero for 
6 = 0, tumbling cannot be initiated unless the initial 
value of @ is finite. However 6 = 0 is never attainable 
in practice for many reasons, and 6 will build up when 
C/A is less than unity. By differentiating (22), 6 can 
be shown to have a maximum at @ = tan (1/+/2) 
= 35°20’. Due to cos’ 6, 6 will diminish to a small 
value near 6 = 90°, and an infinite time will be re- 
quired to reach this angle. 

For small values of 6, (22) is approximately equal to 


C= 160 (23) 


and the time required for the attitude angle to change 
from 6) to 6; 1s 


— 2 
a " ms ; (24) 


Numerical Example 


Let the two solid disks be aluminum, 4 inch thick 
and 24 inches diameter, and the flexible tube be 0.032 
inch stainless steel 6 inches in diameter and 24 inches 
long. The quantities required for the computation of 
K are: 


C = 8.16 Ib in. sec” 

A = 20.44 lb in. sec” 

ne =. 0. 1136 lb in. * see” 

Ve Aen 

£ = 3.0 1n, 

[=F 2/7 ine 

(f= PAD nae 

B29 G10 thin, © 
Assuming y = 0.05 and w = 507, the value of K is 
662 10°. Thus for the body to undergo an attitude 


ete ‘5 a : : 
change from 1° to 10°, the time required, as calculated 
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: 
from (24) is 


t= — X 10° = 3480 sec 
= 58.0 min 
V. Example Including Inertia Forces 


In the preceding section, inertia forces due to th 
elastic motion were neglected, and so no resonan¢ 
phenomenon is observed in the solution. A proble1 
where inertia effects are important is the followin 
(Fig. 5). 

The rotating cylinder of radius F is rigid. Four un 
form dissipative, elastic cantilever beams project frot 
its center at right angles to the axis of symmetry. Th 
beams bend only in the z direction. Assuming the elasti 
deformations w(x, t) to be small compared to the rigi 
body motion, the equation of motion of each beam is 


Bie ae 


Le? (4 


where a, is given by (15) and & has been replaced b 
x + R (Fig. 5). ET is the flexural rigidity of one of th 
beams, JV is the beam mass, L is the beam length. 

The deflection w(x, ¢) may be expressed as a serie 
of products of generalized coordinates q,(t) and th 
normal functions y,(«7) of the beam, 


(2 i es ¥ tal t)Yn(x) (26 


n=l 


The normal functions y,(2) are assumed normalize 
to the beam length, i.e. 


L 
i PRD pia (27 


Substituting (26) into (28) and using the fact that th 
normal functions are orthogonal the equations + 


Fig. 5. Spinning Cylinder with Projecting Beams 


| 
| 
| 
| 
| 
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Fic. 6. Variation in the Rate of Tumbling including 
esonance Effects. 


motion become 
L 
in (t) <2 aE ts 7 =) Oo Qn(t) — a, V(x) Aw dx (28) 
Qa L Jo 


where a complex structural damping factor zy/27 has 
been assumed, and , is the nth natural frequency. 

ter substitution of a, from (15), the integral in (28) 
an be evaluated from the tables in [7] using the con- 
stants a, and B, tabulated in [8], to give 


peelaee 


Pe Aon C 
PALe, : 


(29) 
(1 + Ra, B,) sin dt sin 6 cos 6 


Solving this forced vibration equation and integrating 
he damping force for each mode over one precession 
yele, the over-all dissipation per cycle is given by 


32C’M yw sin’ 6 cos’ 0 os 
A? IDF 0—1 


(an Batt sr 1)’ (30) 


C-O Oye 


where the dissipation has been increased by a factor of 
four because there are four beams. Since the quantity 
3, increases very rapidly with increasing n, all but the 
fi ‘st term of the above series will be neglected. Then 
sing (3), (8), and (16) the equation for the cone 


= 


angle drift rate is 


6 — 
16CM (a, B, R + 1)’ wy sin 6 cos’ 6 ney 
DN 2) 
A? ‘1703 1 -[( — CV er Iy+ 
qT By 1 l mi 9 cos 0 = 
or 
Ae K’ sin 6 cos’ 0 
a —p cos 6)? + Gar Oe 
Qa 
where 


er 


Note that the resonance term in the denominator of 
(32) allows 6 to become large if the forcing frequency ¢ 
coincides with the beam natural frequency , for some 
6. This effect will be present for all p > 1. The effect is 
illustrated in Fig. 6, which is a plot of (32) for the 
case (y/2r) = 0.01, p = 1.2, 1.6, and 2. Note that the 
envelope of the peaks in Fig. 6 is identical in shape to 
Fig. 4. 
The solution to (32) is 


i= E ae Cl (sec 86 — sec 6) 


0 
2 wena = 
os f 4. eo +p — 25 | In eres 


0 
tan — 
2 


+ p'(cos 6 — cos 8) 


Equation (34) is plotted in Fig. 7 for the case (y/27) 
= 0.01, p = 1.2, 1.6 and 2.0, 6(0) = one degree. The 
resonant behavior in these cases causes the angle to 
diverge quite rapidly. For some eases involving lightly 
damped resonance the assumption that @ is slowly 
varying, made in this analysis, may not hold. 


CONE ANGLE @ DEGREES 


K't DIMENSIONLESS TIME PARAMETER 


Fig. 7. Attitude Change of Unstable Body including 
Resonance Effects. 
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VI. Conclusions 


1. The attitude change of an unstable, dissipative 
body can be predicted by calculating the oscillating 
stress field in the body, and relating it to the energy 
dissipated. 

2. The acceleration given by (14), used as inertia 
force per unit mass, and used in the appropriate dif- 
ferential equation for the structure, will result in the 
amplitude of cyclic vibration from which the rate of 
energy dissipation can be determined. 

3. If only the z component of acceleration is involved, 
as in the illustrative examples, the equation for 6 will 
have the form of (32) where K’ will take on different 
numerical values depending on the configuration. If 
inertia forces due to the elastic motion are neglected, 
the equation for 6 will have the form of (22). 

4. The assumption made in the derivation of equa- 
tions other than (23) and (24) is that 6 is small. This 
assumption is justifiable in most cases, but may be 
violated, for some cases, if resonance is important. 


VII. Nomenclature 


A,C = moment of inertia about pitch and spin 
axes respectively 

h angular momentum vector 

6 = spin axis attitude measured from h 

¢ = spin rate 

v 

1 


= precession rate 
2,0 = body axes 


w;j;%, = angular velocity about body axes (1, 2, 3) 

ao = initial angular velocity for @ = 0 

du = kinetic energy 

o = maximum cyclic stress 

E -= modulus of elasticity 

y = hysteretic damping factor = 27 X solid 
damping coefficient 

to = cyclic time 

é = coordinate along axis | 

Z = coordinate along axis 3 

V = volume of stressed material 

a = acceleration vector 

do = acceleration of origin of axes 1, 2, 3 

a’ = acceleration relative to axes 1, 2, 3 

v’ = velocity relative to axes 1, 2, 3 

@, @ = angular velocity and acceleration of co- 


ordinate axes 1, 2, 3 
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Me = gyroscopic momnet on one disk 
my = mass of one disk 


m = 2m, = total mass (of 2 disks) 

l = length of flexed tube 

K = constant defined by Equation (22) 

R = radius of spinning cylinder 

NE = beam cross-section moment of inertia 

M = beam mass : 

L = beam length 

w(2, t) = beam deflection c 

7 = nth natural frequency of beam 

y,(a) = nth normal function of beam, correspondir 
to 

g.(t) = nth generalized coordinate of beam defle 
tion 

i =vV-1 | 

a,, Bn = constants for uniform beam, tabulated | 
Reference [6] 

Ise = constant defined by Equations (31) ar 
(32) 

p = frequency ratio defined by Equation (33 


Dots denote derivatives with respect to time. 
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Abstract 
If the angle of impact is chosen to be perpendicular to the 
local lunar surface, (i.e. normal impact), for purposes of 
implifying terminal guidance of lunar vehicles, then it is 
ssible to show, on the basis of a rocket moving under the 
gravitational field of the Earth and a circularly moving but 
a Moon, that normal impact region is restricted to the 
to 


‘ 


ding semi-circle of the visible disc. Further, it is possible 
to show that higher energy orbits will strike closer to the 
center of the visible lunar disc than lower energy orbits. 


Discussion 


A lunar orbit is characterized by the flight time re- 
quired for the rocket to impact the surface of the Moon. 
Clearly, a flight time of 1} days implies a higher orbital 
energy than a flight time of 3} days or 4 days. If the 
condition of normal impact, that the angle of impact 
of the vehicle be perpendicular to the lunar surface, 
be invoked, then, it will be shown in this paper, higher 
energy orbits strike closer to the center of the lunar 
dise than lower energy orbits. Further we designate 
areas on the Moon to be reached by elliptical, hyper- 
bolic, or parabolic orbits for normal impacts depending 
upon the initial conditions near the Earth. 

To show this, we assume: 

(1) Normal impact on the lunar surface; 

(2) Spherical non-accelerating Earth; 

(3) Massless Moon moving in a circular geocentric 
orbit ; 

(4) Two-body geocentric motion of rocket; 

(5) Orbit plane of the Moon coincident with the 
equatorial plane of the Moon; 

(6) A translating but non-rotating Moon; 

(7) At final impact conditions, the distance from 
rocket to Earth is the same as Moon to the Earth; and 

(8) Rocket burnout at “reasonable” altitudes and 
speed. 

By considering both the geometry and dynamics of 
the Earth-Moon-rocket configuration, together with 
the above-listed assumptions, equations involving the 
lunar impact longitude (for the coplanar case) and 
longitude and latitude (for the three-dimensional) 
san be derived in terms of gravitational constraints. 
Sertain regions on the lunar surface can be designated 
is “elliptical regions,” others as “hyperbolic regions” 
und the boundary between the two regions designated 
is “parabolic region.” In order to find limits for the 
iormal impact latitude, it is necessary to find the 
maximum value of p, the semi-latus rectum, by use of 
Shysical constraints. The results of the analysis will 
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Fig. 1. Relative velocity of moon and rocket 


show that the effect of p upon longitude, serves to 
“fuzz out” the sharp separation of the areas found to 
exist for rectilinear cases, but the general characteris- 
tics still hold; i.e., higher energy orbits impact. closer 
to the center of the lunar dise than lower energy orbits. 

This analysis, taken together with the study of the 
effect of retrograde and direct orbits upon latitude will 
also show that for retrograde orbits the hyperbolic 
region assumes one shape (the remaining area occupied 
by the elliptical region); while for direct orbits the 
hyperbolic region assumes another shape (with the 
remaining area occupied by the elliptical region). 

First, consider a simple model of massless rocket 
impacting normally at a longitude A, upon a massless 
Moon moving on a circular orbit. Both are traveling 
under the gravitational influence of a spherical, non- 
accelerated Earth. The justification for the neglect. of 
the Moon’s gravitational influence upon the rocket is 
that the effect of relative geocentric velocity between 
Moon and rocket is the principal effect, while lunar 
perturbations are of second order. Lunar rotation is 
also neglected so that the angle of impact is calculated 
in a non-rotating coordinate frame. 

From Fig. 1 we see the different roles played by the 
velocity of the Moon and of the incoming rocket. Let 
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the following quantities be defined: 


i, = geocentric velocity of the rocket 

iy = geocentric velocity of the Moon 

* = relative velocity of Moon with respect to rocket 

h = longitude of the impact point measured posi- 
tively eastward from the Earth-Moon line 

@ = latitude of the impact point measured positively 
northward from the Moon’s equator 

rm = distance of rocket to Harth 

y = true anomaly measured from perigee positively 
counterclockwise to the object 

& = geocentric speed of rocket 

$m = geocentric speed of Moon 

§ = relative speed of Moon with respect to rocket 

p = semi-latus rectum of geocentric orbit in Harth 
radi 

a@ = semi-major axis of geocentric orbit in Earth 
radii 


The units adopted here are “normalized’”’ units; 1.e., 
the unit of time is a ke’ minute (13.447 minutes) and 
unit mass is the mass of the Earth, and the unit of 
length is one Earth equatorial radius. 

The requirements that impact be normal to the 
surface of the Moon implies that the velocity vector 7 
must pass through the center of the Moon. Also, since 
the Moon’s rotation is neglected, we note that every 
point on the surface of the Moon partakes of the same 
velocity as the velocity of the center. 

Figure 2 illustrates the assumed rocket trajectory. 

By treating the Moon as a sphere and not a dise in 
order to study latitude of impact, the following situa- 
tion is depicted by Fig. 3 where an additional defini- 
tion is introduced, namely J, the angle between the 
geocentric orbital plane and the Moon’s equatorial 
plane. 

In Fig. 8, the only vector that lies in the orbital 
plane of the vehicle is 7 . Since the rocket’s orbit plane 
must contain the center of the Earth, the radial ecom- 
ponent of 7 is nearly along 7. The transverse component 
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Fig. 2. Assumed trajectory of rocket 
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Moon's. 
Motion 


Fic. 3. Angle between geocentric orbital plane and moon 
equatorial plane. 


? 
of 7% however lies in the vehicle’s orbital plane, paralle 
to the k, J plane, and can be resolved into two furthe 
components, one along J and another along k. Th 
normaley condition requires that 7 pass through th 
geometrical center of the Moon. The Moon’s motion 1 
a circular orbit about the Earth is along the negativ 
J axis. 

Three vector equations represent the velocity of th 
Moon with respect to the vehicle, and of the Moon an 
vehicle with respect to the Earth as follows: 


zr = icos¢cosX\ + Jscos¢@sinrA + kssing (la 
Ty = —J&u (1b 
iy = —tiy — J(r)s cos I — k(ré)g sin [ (le 


For the requirement that * = 74 — 7%, , the correspond 
ing scalar equations representing the motion of th 
Moon with relation to the rocket will then follow: 


§ cos @ cos A = fy (2a 
§ cos @ sin dA = (7), cos I — Sy (2b 
§sing = + (rd), sm Jt (2¢ 


In terms of the geocentric orbital constants, the es 
pressions for ¢ and } are as follows: 


Weak N= | st (4/2 cos r) — 1]/ 
tang = + (4/ : sin be 


An upper limit of p, the semi-latus rectum, exist 
from purely practical restrictions. The value of pe 
determined wniquely from initial burnout condition: 
In fact, in the normalized units, p is equal to the squal 
of the angular momentum. The amount. of angulé 
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Earth 
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Fic. 4. Projected angle of rocket from local horizon 


momentum imparted to a vehicle initially is limited so 
that p is certainly limited. This may best be seen by 
‘the following considerations. From Fig. 4 it can be seen 
that: 


y = initial angle of projection or path angle of the 
rocket from the local horizon 
é’ = time rate of change of the trwe anomaly 
Yz.o. = distance of rocket at burnout from the center 
of the Earth. 


(rd) z.0. = Sz.o. COS y. 


But by conservation of angular momentum we have: 


(po = pe constant 
Pp = T3.08B.0. cos? 
2 .2 
Pmax = (Te-0. ax $5 6!) ae 


For physically realizable missile systems, we might 
put as upper bounds for rz.0, and & the following: 


(73.0.)max = 1.3 Earth radii (i.e., 12U0 miles above sur- 
face of Earth) 

= 1.5 radii/k,’ min (i.e., 39,000 ft/sec) 

(1.3)? (1.5)? < 4 Earth radii. 


Tmax 

Pmax 
It is noted that present-day rockets do not burn out 
at altitudes of 1,200 miles vith speeds in excess of 
39,000 feet/sec. Thus, for reasonable lunar trajectories, 
p is less than 4 Earth radii. P/r is then less than 0.1. 
Since 


I 


and, 


ee Pmax , 1 
(ria = / - VE 
(rd) a < Su avaek 


It can be seen that sy is always greater than (7é)a cos J. 
This fact is useful in removing the sign ambiguity in 
equations 3a and 3b by the following argument. 


Let the angle of inclination, J, of the relative orbit 
planes always be measured at the ascending node de- 
termined by a direct orbit (inclination is less than 90°). 
All inclination angles hereinafter are always measured 
at this node so that J may take on values from 0° to 
360°. (The normal convention in astronomy is to meas- 
ure IJ always at the ascending node.) Consider that in 
Eq 2b 


(A) & is always positive; 
(B) Since —90 < ¢ < +90, cos ¢ is always positive; 
(C) &% > (rd)a cos I. 


Thus, the right-hand side of Eq. 2b is always nega- 
tive. Neglecting the case of intersection of the Moon’s 
orbit after apogee, so that 74 is positive, we can con- 
clude that \ is restricted to values between —90° and 
0°. In Eq. 3a, the plus sign is valid. 

If Eq. 2b is divided by Eq. 2c, one obtains 


(rv), cos I — Sy 
(rd)a sin I sin X 


cot ¢ = 


The numerator is always negative, as shown above. 
Therefore, 


(negative quantity) 


cot ¢ = = fe F : : 
g sin I (positive quantity) (negative quantity) 


(positive quantity) 
sin [ 


The extreme value of p = 4 determined previously is 
one consideration needed to calculate the maximum 
value of ¢ obtainable by a ballistic missile. Other values 
of ¢ and X, for combinations of p, a, and J, are given in 
Table I, and Fig. 5A and 5B show graphically the 
possible variation of @ and d with a. 

Note that in Fig. 5A, since physically realizable 
rocket configuration limit p to 4 radii, under the most 
favorable circumstances the maximum absolute lati- 
tude obtainable is on the order of +15° for normal 
impact. In Fig. 5B, the upper branches of the family of 
curves do not appear because the angular momentum 
of the Moon is always greater than the angular momen- 
tum of the rocket. Also note that by various choices of 
p between 0 and 4, and J, a certain amount of overlap 
of normal impact areas is possible. 

Figures 6, 7, 8 and 9, obtained by plotting the values 
given in Table I, show the areas where normal impacts 
may occur for various trajectory characteristics and 
also show the locus of all normal impact points for the 
parabola. The area is over a restricted portion of the 
Moon’s surface, bounded by longitudes of —90° to 0° 
and latitudes of about +15°. The areas for longitudes 
from —90° to —35.3° can be considered effectively as 
elliptical (with the exception of the inclination effect 
which serves to shift the areas to the right of —35.3° 
for inclinations between 270° and 360° and 0° to 90° 
and to the left of —35.3° for inclinations between 90° 


The Journal of the Astronautical Sciences 37 


| 
faze 
+10 | I = 90° 
ee See l ell Se 
| I = 90° 
ee ae ee oe ee ee 
PE ersten 
p=0 | re 
Race eS if r a=5 a=r a 
p=l = 270° 
m7 5 [| 
oO 
p=4 I= 270 | rae 
ja 
| oe ioe 


- 15 


| 


Fra. 5B. Plot of longitude vs. semi-major axis 


and 270°), and the areas from —35.3° to 0° as effec- 
tively hyperbolic. Also, as can be seen from Fig. 6 
through 9, for what is usually termed as direct orbits; 


1eenO eae P< O08 2270 = P< 360° the iw perbolic 
regions take one shape while for what is conventionally 
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TABLE I 


p = 0 (Rectilinear Orbit) 


a Tm +rM/2 +0 —IM 0 
|  —4520 —90°0 —35°3 =30°0) |) 3 
p=1 
7 
Tu 
a TM ae too ™ 0 ; 
Iu j 
I r ? x re) nN fe) en ¢ r q 
0-|—4122] 0° |—90°) 0° |—3228] 0° |—26°8| 0° | 0° 
90°|—45°0| +5°1/—90°| +7°2|-3593)+-4°2| 3020 +3°6) 0° 0 
180°|—48°4| 0° |—90°| 0° |—38°5| 0° |—33°0| 0° | 0° O 
270°|—45°0| —5°1/—90°| —7°2|-35°3| —4°2/-30°0|—3°6| 0°| O 


I » ra) » ’ » ? r te) A | @ 


4 0°/—36°0} 0° |—90°} 0° |—27°9} O° |—23°4) 0° | 0° @ 
90°| —45°0| + 10°2|—90°|+14°2| —35°3) +-8°2|—30°0|/+722) 0° O° 
180°|—51°3} 0° |—90°| 0° |—41°5| O° |—35°8) O° | 0° O 
270°| —45°0} —10°2| —90°| — 1422) —35°3| —8°2) —30°0|—7°2) 0°| O° 


termed as retrograde orbit; i.e., 90° < I < 270°, the 
hyperbolic region assumes another shape. p = O is not 
shown, since on these figures, the plot will merely be 4 
straight line falling upon the equator. 
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Error Analysis of Satellite Orbits 
in the Presence of Drag 


a i at eae e 


Frederick V. Pohle* 


Abstract 


First order perturbation solutions of satellite orbits in the 
presence of drag were used to determine the error in the 
altitude, at the end of a half revolution, due separately to 
(i) errors in the initial altitude, (ii) initial horizontal speed, 
(iii) perturbation parameter, and (iv) decay parameter in the 
exponential density function. The error in the speed, at the 
end of a half revolution, was also determined as a function of 
the same parameters as in (i) to (iv). The present analysis is 
restricted to the case of a horizontal velocity at the initial 
position. In the absence of drag, information can be obtained 
easily on the effect of initial non-zero vertical speed upon the 
orbit. 


I: Introduction 


The trajectory analysis of a manned satellite may 
be divided into three stages: (1) take-off, (41) orbital 
phase, and (iii) re-entry. This paper is concerned only 
with (ii), namely, it is assumed that certain conditions 
are prescribed at the end of (i) which serve as the 
initia] conditions for (ii) and that these are sufficient 
to allow the satellite to stay in orbit for some time 
before the re-entry phase occurs. 

Great accuracy in guidance is required to place the 
satellite in a given orbit. Inevitably there will be devia- 
tions (errors) from the assumed initial conditions. It 
is important to be able to determine the effects of these 
errors upon quantities of physical interest. In particu- 
lar, the initial altitude, initial speed, perturbation 
parameter (which involves mass, drag coefficient and 
area of the satellite), and the decay constant in the 
assumed exponential density, may all vary from as- 
sumed or given values. This paper is concerned with the 
determination of the effect of these errors upon the 
altitude and speed at the end of a half revolution. 

A precise orbital analysis would require the con- 
sideration of many effects. However, in phase (ii) the 
effect of drag is most pronounced and simplifying as- 
sumptions can be made; the precise assumptions are 
stated under Equations of Motion. A more complete 
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discussion is contained in [3] and the re-entry phase i 
discussed, for example, in [2]. The present work is at 
extension of the work of Roberson, [1] who obtainec 
first-order drag corrections. 

Roberson, [1], was concerned only with the sotellil 
life-time study. Rather than compute each circut 
numerically on a digital computer, each complet 
circuit was computed analytically and the termina 
conditions of one circuit were used as the initial condi 
tions for the next circuit. In effect the orbital analysi 
was based upon difference equations in the revolutiol 
number and only the difference equations were evalu 
ated on the computer. This approach to the study o 
satellite life-times is a significant improvement ove 
purely numerical studies. In the present problem it 1 
of more importance to analyze each half-revolutiot 
near re-entry and this has been done. 

For completeness, the method of solution is out 
lined and in particular all quantities which occur in the 
final calculations are fully defined. The error deriva 
tives have been calculated from the first-order drag 
solutions under the uniform assumption of a zero initia 
vertical component of velocity. Since an error analysi: 
for the vertical speed is also of value, a simplifiec 
analysis is presented for this case. The present analysi: 
must be modified if a drag correction is also to be ob 
tained in this case, and is not attempted in this report 

Finally, some additional remarks have been made 
on the possible presence of lift, second-order drag ef 
fects, non-zero initial vertical speed, and other modifi 
cations in the analysis. 


II: Equations of Motion 


The earth is assumed to be a sphere and the atmos 
phere is assumed to be non-rotating. Under these con 
ditions the orbit will be in a plane. If the drag is tan 
gential and proportional to the local density and to th 
square of the speed then the equations of motion i 
polar coordinates (7, B) are: 


dr/dt’ — r(dg/dt)? + K/r’ 
= —(CrA/2m)p(h)r[(7)” + (78) 7? & 
r(d’B/dt') + 2(dr/dt)(dB/dt) 


—(CpA/2m) p(h)rB[(r)* + (7B)? 2 


It will be convenient to introduce the appropriat 


n-dimensional coordinates as 
f= a/r (3) 
n = Ka/(r°B)° (4) 


The Equations (1, 2) contain r(¢) and 6(t): through 
3, 4) the variables will be changed to (8) and 7(8) 
ud primes (’) will be used to denote differentiation 
vith respect to @. If the substitutions (3, 4) are carried 
mut in (1, 2) the resulting equations, following Rober- 
son [1], can be written as 


. 


Zz gE = y (5) 
3 a'/n = (90/8) + (8/8)? (6) 
Equations (5, 6) are the equations used in the subse- 
juent analysis. 

111: Perturbation Solution 


~ It is assumed that expansions can be written in the 
‘orm 


l| 


E(B) = E° + (vp)t? + --- (7) 
n(B) = 9° + (vpo)n? + --- (8) 


where the non-dimensional parameter (vpo) is assumed 
(0 be small compared with unity; only first order cor- 
ections will be computed and the zero-order terms are 
zaturally the Kepler solutions which define the refer- 
nce (drag-free) orbit. The zero-order equations are: 


ee ae By _ Ay! (9) 
ee = (10) 


ind it follows that 7° is a constant, which becomes the 
ight side of (9). 

As initial conditions it will be assumed that at ¢ = 0, 
= 0 and that r = 1(& = &); V, = O(2 = 0) and 
1 = (Ka/r;V,') where Vz is the given speed. The 
olutions to equations (9, 10) can then be written as 


1 (8) = m = constant (11) 
E(B) = ml + a cos 6] (12) 

here 
a = (f/m) — 1 (13) 


he eccentricity of the orbit is given by (13), and 
11) expresses the fact that the area is swept out 
niformly in time. This Kepler solution is used in the 
rst-order drag solution. The equations for the cor- 
ection are: 


ae a ay = 1” (14a) 
(0) (0) (0)/ 2771/2 
Fe (3) p(é”) £'(g) 
a os E i (En ) | fue 


and 


1) 
qo (lion) i plm(1 + 4 cos B)] x 
(15) 


+ 2a cos 6 + 2 dB 
(1 + «cos 6)? 
Equation (15) follows from (14b) and (12). It is 


necessary to expand the terms which do not involve 
the density, as follows: 


1 9 21/2 oO 
\ “3 om oe = 2 ¢n(a) cos nB (16) 


The coefficients ¢,(e) are defined in the Notation and 
a sufficient number of terms are shown to make possible 
solutions which systematically contain terms up to the 
third order in the eccentricity. It should be pointed 
out that the third coefficient of go(e1) in [1] is given in- 
correctly. 

The density function in (15) can be written as 


p[lm(1 + « cos B)| = p(m) exp (A cos 8) (17) 


where 
A= mare (18) 


and p(m) is defined under Notation; clearly, at 6 
= (4) the exponential term in (17) is unity and the 
density p(m) is that at the altitude mid-way (in angu- 
lar measure) between perigee and apogee. It is stated 
in [1] that the value of \ is about 25 for qg = 0.2. A 
value of the eccentricity of s49 would correspond to an 
altitude change of about 400 miles; in this case \ would 
be about 6.25. The only density which occurs in the 
final calculations is p(m); the actual densities vary 
from this value by the factors e* and e°. 
From [4], page 164, No. 103, it is known that 


i exp(\ cos 8)cos nB dB 
Tv 
; (19) 


== a exp(A cos 8) cos nB dB = I,,(\) 
T /0 


The result- (19) makes it possible to write the results 
for the half-revolution as well as for the full-revolution 
and attention will be restricted to the former case. 
The result for (15) can be written as (8 = 7) 


a (ar) = [ep(m)/pol Dd only (20) 


where the argument of ¢, is « and that of J, is X. 

Since the right hand side of (14a) is known the 
solution can be written immediately in the integral 
form 


i> 2) 


£(8) = [o(m)/pol >) en(a) X 


m—0 


B 
i exp (\ cos B) cos ny [cos(B — y) — 1] dy 
0 
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and therefore 


E(r) = —[mo(m)/2po] oe esa ia ee) 


or 
EY Cr) = [mp(m)/po] 2, enBn (23) 

where 

Isp = B,(X) = = (Ce =I aL; = 1 Pe and 


B, = B_n 


The first order corrections are given by (20) and 
(23). In the following it will be necessary to evaluate 
I,’ and B,’. From [4] (page 163, No. 90), it is known 
that 


Ln = (5) (Let Ingt) (24) 


If the definition of B, is used, together with (24) 
the result is that 


B,! = (3) (Bra + Bays) (25) 


which is the same functional form as (24). Thus (20) 
and (23) have the same form and the results for one 
will be the same as for the other in terms of B, and 
I, . The differentiations are taken with respect to the 
variable i. 

It is now possible to write down the results for € and 
n at 8 = 7, as follows: 


E(1) = (2m — £1) =F mvp(m) D1 mB, (26) 


n(n) = m+ rm) Miele (27) 
where, for convenience, it is to be recalled that 
& = (a/n) 
m= (ga/Ve)& = fb 
$1 = (ga/Ve) (28) 
a = [(&/m) — 1]; on = en(a) 
A = mare ; ‘Ee PONG Bie Ban) 


The first-order solutions (26, 27) will be used to- 
gether with the definitions listed in (28) to determine 
the error derivatives. 

IV: Error Derivatives 


The quantities which will be varied are & , ¢, a1, v. 
If the first alone is varied, for example, this implies an 
error in the initial altitude. The derivatives which are 
needed in the final calculations are 


dm/0& = 2m/k 
da /0& = —1/n 
0A/d& = ay[l — (2m/é)] 


dp(m)/d& = (Qarm/é) p(m) 
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dm / 96 ie £” 


da/0t1 = —b'/m ( 
Od/dt, = anky 
dp(m)/db. = —anki p(m) 

(8f1/9Vs) = —29a/Ve° (3 
dp(m)/day = meh) (3 
d\/doy = mea = (& — mr) 7 

} 


The rate of change of final altitude with respect t 
the initial altitude can be written as 4 


- 
OE(r) _ me 
aaa [4(m/é&) — 1] 
+ mvp(m) 2 ben Bn! call — 2(m/Ei)| : 
i 
+ ¢,'B,(—1/m) + mvp(m) on(2m/é) a ¢n B 
where 
gn’ = do,n/de = Pn(a) 
iByee — dB,,/dX = (3) [Bea oe Brasil 
dE(r) _ Pa 
cn [4(m/&) — 1] 


ine} 


aia mvp(m) ys , B,(—1/m) an (2a m/é) Pn BS (33 


v= 


+ = [1 ar 2(m/&)]¢en( Bra —- Baan 


2 


ee = a + varp(m) a [en B,!(—a &’) 


sp Pn Dn ( ar £°/m') | 


+ (vr) p(m) ay Ey ey (Gn Ba 


n=0 

ag ‘ 2 (34 

(a) = 28° + vep(m) do 

Oo71 n=0 

‘- Eb /m ®,, Be a ay Ey Pn iy 
= (ay EPS) Bee se Baws) 

a& (mr 2 
oe = mvp(m) Dy en Br'(& — m) 

Oa n=0 


+ mvmp(m) dD) on Bn 
n=0 (35 


= vrp(m) ay: {(& = m) (1/2) 


n=( 


OrABra ap Bratl aI ™m Pn B,} 


I 


mp(m) Pe nBn 


(Cp Aa/m) 


nd changes in the values Cp , A and M can be related 
0 changes in ». The derivatives which show the effect 
f different quantities upon the altitude are given by 
33) to (36). 

The corresponding derivatives for the speed can be 
written down directly when it is recalled that (26) and 
27) are essentially the same in the drag terms. 

Let 


(m7) = n()/E(r) 
L/§ (a) JOS (ar) /OE, 
= [1/n(m)]on(x)/d& — [2/E(a)]0E(r)/0& (37) 
1¢()|9¢() /dh1 
= [1/n(m)]On(m) /d6. — [2/&(m)]0E(r) /OK 
nd 
/6()J9¢ (4) /don 
= [1/n(m)]on(m)/d0q — [2/E(m)J0E(r) /Oar 
[6 (a) \0¢ (4) /av 
= [1/n()]on(r)/dv — [2/&(m)]0E(m) /dv 


The definition (37) can be used to determine the 
erivatives which have been written in (38). The 
erivatives in the second terms on the right hand sides 
f (38) are known from equations (33) to (36). The 
erivatives in the first terms can be obtained from an 
<amination of (27) which shows that the drag terms 
ill be the same as those of (26) if B, is replaced by 
,; therefore the only term which is left is the first 
rm of (27) which is the constant term from the 
ference orbit solution. The derivatives of the first 
rms are known from the first results of (29) and 
30), respectively, and therefore the derivatives for 
1e speed as a function of initial altitude, initial speed, 
rag parameter, and exponential decay constant are 
so known. 

Examination of the leading term of (33) shows that, 
the case of a circular orbit, the drag-free term is 3, 
at is, an error in initial altitude is multiplied by 3 
the end of a half revolution. This can also be proved 
rectly and it is easy to understand the result: when 
is recalled that the speed which is appropriate to a 
ven altitude for a circular orbit will cause the satel- 
e to move in an ellipse of small eccentricity if the 
eed is kept constant as the altitude is changed 
ghtly. 

The present analysis must be substantially modified 
initial vertical speed is to be added as a parameter. 
1e presence of such a parameter causes the apse line 
the reference ellipse to be rotated and to introduce 
ditional terms which will not modify the analysis if 
‘ull revolution is considered but which will change the 


(36) 


= 
I 


(38) 


analysis completely if only a half-revolution is con- 
sidered. If the analysis is based only on a drag-free 
orbit and if the initial speed is given through the value 
&’ then the equation of the orbit is 


E(B) = m + [(& — m)] cos 8B + &’ sin B 
and the eccentricity of the orbit is given by 


({(&/m) = ie + Guenae” 


If, for example, the vertical speed was initially zero, 
so that the orbit would be circular, an error in this 
value would cause the orbit to be slightly elliptical and 
the errors in altitude, for example, from this case 
alone could become appreciable. 


V: Concluding Remarks 


The present analysis has neglected certain effects 
which may also be important. For example, the presence 
of lift in the governing equations would seriously modify 
the results; this additional force can be included but 
the analysis must be modified. If the lift-to-drag ratio 
is small then a reasonable approximate solution is 
again possible. 

In the work of [1] the density is assumed to be small] 
enough for the quadratic approximation to be negli- 
gibly small. In the present case where the flight is at 
lower altitudes the second-order correction could be- 
come significant. The higher terms can be included 
systematically in principle but the practical difficulties 
are such that a modified analysis would have to be 
used. 

Finally, other basic equations could be used with 
profit in the analysis. In particular other dynamical 
variables can be used and different perturbation param- 
eters can be introduced to fit special circumstances. 
Energy relationships can also play a significant role in 
the determination of the variation of orbital elements. 


Notation and Definitions 


a = radius of the spherical earth (ft) 

A = cross-sectional area of the body, normal to 
the flow field and assumed constant (ft”) 

Ba >= =(4) Una SF 20 + Teas 
BL = BG 20) B, = BO 

B,! = (ey CBeet aR Bry) (n = 0)B,’ 7 dB,,/dx 

Cp = drag coefficient (constant) 

g = acceleration due to gravity at surface of the 
earth (ft/sec”) 

h = altitude above the surface of the earth 
=r—o (ft) 

ih = [,(X) modified Bessel function of the first 
kind: } = amea;l. = itn 2 0) 

Pe oe (2) 4 La) (0 0) Ty vert, aX 

Kk = ga (f{t°/see’) 

m = mass (slugs) 

n =-intever, 0, V2, °-- 

ip = distance from the center of the earth to the 
orbit (ft) 
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ry 


Ve 


by 
cT 


44, 


initial value of r (¢ = 0) (ft) 

time (sec ) 

velocity component perpendicular to the 
radius vector r (positive in the sense of 
increasing £) (ft/sec) 

velocity component parallel to the radius 
vector r (positive outward) (ft/sec) 

d/dt; differentiation with respect to the time ¢ 

altitude parameter for the reference orbit 

angular advance of the body from a fixed 
reference line (inertial axis). The location 
8 = 0 corresponds either to perigee or to 
apogee, that is, it is at an apse. (radians) 

constant in the exponential density function 

eccentricity of the reference (drag-free) or- 
bit = (aym) 

a/r; non-dimensional measure of the altitude. 
The quantity is always less than unity and 
a change of one part in 4000 is equal to a 
change in altitude of one mile. 

a/r; initial value of &. 

drag-free value of &, as a function of 8 

first-order drag correction to &, that is, 
E = £(8) = E? 4 (vpré” 

(ga/Vae)é 

initial value of y 

drag-free value of 7 

first-order drag correction to 7, that is, 
n = 9(8) = 9° + vpoq” 

ga/Vg non-dimensional measure of V4, 

initial value of ¢ 

density at the altitude h. Typical values are: 


(slugs/ft*) 
h(feet ) (slugs ft”) 
300,000 6x 10° 
400,000 6.5 x10 
500,000: 4385x1057 
600,000 45 x 10% 
700,000 (3 erion 
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The atmosphere is assumed as locally 
ponential, that is p = exp (a& + 6) wh 
a and 6 are known constants and exp ( 


=. 
po = density at the initial altitude (71 — a) 
(slugs/ft 
o(m) = reference density at the altitude correspon 
ing to 7, that is the altitude = (a/m) — . 
This altitude is intermediate in value 1 
the maximum and minimum altitudes. ; 
(slugs/ft 
go =1+4 (8)(a)’ + (4)(a) + --- i 
a = —(a) — ($)(a)’ — (&) (a)? + °° 
ge = (4)(a)? + Ge) (a) + 
g = —()(a)t-- ; 
O4 = —(44)(a)' + men 
by = (dgo/da) = ($)(a) + 48)(a) +> 
6 = (dg/da) = —1 — (%)(a)” — (3)(4 
b = (dyx/da) = —(4)(a) + (#)(a) + °° 
@, = (dg/da) = —($)(a) + °-: 
Py = (dy/da) = —(44)(a)° + <-> 
y = (pAa/m; drag parameter (ft*/slugs 
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The Intersection of Coplanar, 
~ Confocal Conic Sections 


Geza 8. Gedeon* 
bstract 


A simple method is presented to calculate the inter- 
sections of two coplanar confocal orbits. 


iscussion 


Given two coplanar position and velocity vectors the 
iestion arises: where will their orbit intersect? In the follow- 


< - 
ea ae os 


Fig. 1. Notation 


g a simple method is presented to calculate the intersections 
these orbits. Using the notations defined on Fig. 1 and in 
e Nomenclature the semi-latus rectums and eccentricities 
the two orbits can be calculated from the following equa- 
ms: 

rK? cos? ¢ (1) 


p 


09 9g a) 


e= V1—p/a (3) 


r the derivation of these equations see [1]. The true anom- 
es related to the given radii are obtained from the equation 


Technical Notes 


line of apsides is 
€ = 0 — (6 —@) (5) 


If ‘& new coordinate is chosen to bisect the angle e, then the 
common radii of intersection can be expressed as 


Po v3 i) 
1 +4 cos (¥ + $) 1 +e cos ( -§) ©) 


Using simple trigonometric relations Eq. 6 can be rewritten as 


[E-S Joos cosy +[£+2]sin$siny = 2 Z (7) 
Pp Pp 2 p Pp 2 po. oP 


If an angle 6 is introduced by the definition 


€ 
t = 
jm tan fee (8) 
3 
where 
ae} 
€0/e 
tae emai 9 
oe (9) 
€0/€ 
then Eq. 7 can be expressed as 
£ cos 6 
~y=6+cos7! €0 € (10) 
cos = 
2 
with 
as Po/p — 1 
pop (11) 


€0/e 


Figure 2 presents the intersection factors € and ¢ as functions 
of the latus rectum and eccentricity ratios. 

In Eq. 8 the are tangent has two values. If these values 
are measured from the polar axis of the chosen coordinate 
system, a full line is obtained and Eq. 10 yields two sym- 
metrical solutions with respect to this 6 axis. 

The true anomalies of the intersection related to the line 
of apsides of the intersecting conic sections can be obtained 
by adding resp. subtracting ¢€/2 from y. 

Eq 10 can be used to determine the conditions when the 
conics intersect, are cotangent, or one is wholly contained 
in the other. For this purpose Eq 10 will be rewritten by Eq 8 
in the following form: 


Ref. [2] 


6=o+ tant | 


Kk? - 1 


om Fig. 1 it can be seen that the skew angle between the 


*Head, Astrodynamics Laboratory, Astro Systems and 


fo) ane e/2 
sec? (y — 6) = = cos? S| ae wre). (12) 
| (4) ; s 
or 
- éo7 | 1 1\1-+ cose 
wr -0-8[L4 (0-2). a 


search Laboratories, Norair, a Division of Northrop 
rporation, Hawthorne, California. 


If € and ¢ are resubstituted then Eq 13, after a little rear- 
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Fra. 2. Intersection factors — and ¢ as functions of latus 
rectum ratio and eccentricity ratio. 


rangement, becomes 


sec? (W — 6) = 


Since the secant must be greater than one, the numerator 
must be greater than the denominator for intersection and 
equal for cotangency as it has been already shown by Lawden 
[3] through geometrical argument. To give a physical ex- 
planation to Eq 14 the following steps will be taken 


tan? (Y — 6) = sec? (YW — 6) — 1 


€¢ 2 | 
CDE ty tee 
P Po ppo_| (15) 


= ae es 

Pp PPo Po 
Both the numerator and the denominator will be then 

multiplied by aappo yielding the following 
tan? (W — 6) 

i —[ao?(1 — eo?) + a2(1 — e?) + 2(ae) (do 9) COS € — 2aay| (16) 
A seb Bs 

PPo 
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If the semi major axis is defined as positive for ellipse 
negative for hyperbola then the linear eccentricity can 
introduced with the same sign convention. 

Thus Eq 16 becomes: 


1 Sore SPOR See 
coz — 2coc cos e + c? — (a — a)” ( 


tan (Wy — 6) = 
‘ (po — py? 
hh = 
PoP 
In this expression 
coe — 2eccose + 2 = fe ( 
f is the center spread. ‘ 


To have real solution for Eq 17 the center spread must 
greater or equal to the semi-major axis difference in case 
similar conics, and must be less than or equal to the sum_ 
the absolute values of the two semi-major axes if one of #l 
conics is a hyperbola. In case of equality y = 6 is the pou 
of cotangency. : 


| 
Nomenclature ; 
wh = G(m + m) & gk? planetary constant : 
r = radius vector : 
Vo = Vu/r = circular orbit velocity : 
V = orbital velocity 
K = V/Vo = Kepler Number 
a = semi-major axis 
c = linear eccentricity 
e€ = c/a = numerical eccentricity 
p = (1 — @&a = semi-latus rectum 
f = center spread 
@ = angle of elevation 
a@ = angle between two position vectors 
6 = true anomaly measured from pericenter 
€ = skew angle between the line of apsides 
wy = true anomaly measured from the bisector of € 
6 = axis of intersection symmetry measured from the } 
sector of € 
polP — 
a ES = first intersector factor 
Po/P 
SEE 5 
€o/e 
Cxoei eee = second intersector factor 
DolD _ 
€0/€ 
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Two Maneuver Ascents to 
Circular Orbits 


J. F. Wolfe* and D. DeBra* 


Characteristic velocity is a common criterion for evaluating 
cent and transfer paths [1]. It will be used in this note 
ough other criteria such as mass ratio may give different 
id sometimes more significant results [2]. 

Large circular orbits may be attained by coasting on a 
ansfer orbit after initial thrust and then performing a 
cond circularizing maneuver when the proper altitude is 
ached. The total propulsion requirements can be estimated 
r evaluating the initial velocity at the start of the coast 
bit and adding the velocity impulse necessary to accom- 
ish the second maneuver. The sum is called v7 , the charac- 
ristic velocity. Use of additional impulses may reduce vz but 
also reduces reliability and requires additional staging for 
ckets not capable of multiple restarts. Therefore, only a 
igle velocity impulse in addition to the initial velocity will 
- considered. 

The second maneuver will occur at apoapsis of the transfer 
bit, which must coincide with a point on the desired circu- 
r orbit. The usual methods for finding extrema show that 
is condition minimizes v7 with respect to the co-flight path 
ile of the transfer orbit at the point of second maneuver. 
his condition for minimization is true for non-coplanar as 
sll as coplanar transfer orbits. 

In this discussion quantities are normalized to the initial 
ius and circular speed at the initial radius, which auto- 


APOAPSIS OF TRANSFER 
ORBIT 


FINAL ORBIT 


TRANSFER 
ORBIT 


Fic. 1. Definition of notation for ascent trajectory to a 
cular orbit. 
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Fic. 2. Velocity requirements for circular orbits 7 = 0 


matically normalizes the constant of an assumed inverse- 
square-law, central-force field. Notation is defined in Fig. 1. 

The characteristic velocity for coplanar ascents is just a 
function of the radius of the desired circular orbit and one 
other parameter which defines what path was followed to get 
there, e.g., 8, 8 etc. The final velocity maneuvers (v. — vq) 
and vr are plotted in Fig. 2. It may seem surprising that v7 is 
not a monotonically increasing function of the final orbit 
radius, 7. At small r the behavior is as expected, but for large 
r the final maneuver decreases more rapidly than the initial 
velocity increases thus causing a maximum in characteristic 
velocity as noted in [1]. 

The expressions for these quantities can be found para- 
metrically in the co-flight path angle by using the conserva- 
tion of momentum and energy along the transfer orbit. 


Vo = 0, 81n B/r 
2(1 — 1/r)/(. — sin 26/r?) 
Circular speed is (1/r)!’? so the characteristic velocity is 
Vp = V4 + (Ve — Va) = Ve + O41 — Vq/0.) 
= (l/r)? + [2(1 — 1/r)d — sin B/r)/(1 + sin gre”? 


It is sometimes more convenient to express the results as 
a function of transfer angle 6. Since the angle @ is the supple- 
ment of the true anomaly of the initial point on the transfer 
orbit, the relation between 6 and @ is 


ctn 8 = esin 6/(1 — e€cos @) (2) 


ll 


9 
Oy 


By solving the orbit equation for the semi-latus rectum and 
by equating the expression for this parameter at apoapsis to 
that at the initial point, one gets 


e = (r — 1)/(r — cos @) 
and from (2) 
sin? B/r? = (1 — cos 6)/[2r? — (2r — 1)(1 + cos 8)] 


which can be substituted in (1) to give the desired result. 
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VELOCITY REQUIREMENTS FOR 
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Fic. 3. Velocity requirements for circular orbits @ = 90° 


OPTIMUM TRANSFER ANGLE 
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Fic. 4. Optimum transfer angle for non-coplanar transfer 


When the desired orbit plane does not contain the initial 
point, the geometry is specified by r and the minimum in- 
clination possible between the transfer orbit and the final 
orbit, tmin. Again injection will occur at apoapsis of the 
transfer orbit. The second maneuver must also change the 
planes, so it must be the vector difference rather than the 
scalar difference plotted in Fig. 2. The characteristic velocity 
is then 


vp = 0; +| 5, —0g| = 05 + (V2 — Qvva cost + v42)"? 


where 7 is the inclination of the transfer orbit to the final 
orbit. 

Choosing the transfer orbit 6 = 90° results in the least 
inclination. This case is plotted in Fig. 3. For some 7, for 
example 30°, extrema of both types occur. 
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MINIMUM TWO IMPULSE VELOCITY 
REQUIREMENTS FOR CIRCULAR ORBIT 


NORMALIZED MINIMUM VELOCITY (V_) 


| 2 4 10 20 40 100 
NORMALIZED RADIUS (r) = 


Fig. 5. Minimum two impulse velocity requirements fo 
circular orbits. 


The minimum v7 does not occur for minimum inclination 
By including the geometrical relation 


siN dmin = Sin? sin 6 . 


Vr may be minimized. 

Then for given r and imin , one additional parameter e.g. 
1, 6, 8, ete. specifies the transfer orbit. The turning poin 
may be computed by setting the partial derivative of vr wit! 
respect to this parameter equal to zero. Here av7/00 = 0 wa 
used to compute 45+ , the value of 6 which minimizes v7 . 

The derivative is-algebraically complicated and not con 
venient to work with except forr = landr = o. For tmin 7 
0 the value of @opt is 90° for r = 1 and as r — ~= it is give 
by the expression 


sin? Oo5pt = 2SIN tmin — SiN? tmin 


Other values of Oop are plotted as a function of r in Fig. 4 
The corresponding values of v7 are given in Fig. 5. 

In the non-coplanar case the optimum value of 6 wi 
generally be less than 90°; that is, the flight path angle wi 
be positive. The requirement for a large flight path ang! 
may reduce propulsion efficiency during initial thrustiny 
Relaxing the constraint of injection at apoapsis will increas 
vr but in permitting a smaller initial flight path angle a 
overall increase in efficiency may result. 
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